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Electrostatic correlation effects in inhomogeneous symmetric electrolytes are investigated within a 
previously developed electrostatic self-consistent (SC) theory (R.R. Netz and H. Orland, Eur. Phys. 
J. E 11, 301 (2003)). To this aim, we introduce two computational approaches that allow to solve 
the SC equations beyond the loop expansion. The first method is based on a perturbative Green's 
function technique, and the second one is an extension of a previously introduced semiclassical 
approximation for single dielectric interfaces to the case of slit nanopores. Both approaches can 
handle the case of dielectrically discontinuous boundaries where the one-loop theory is known to 
fail. By comparing the theoretical results obtained from these schemes with the results of the MC 
simulations that we ran for ions at neutral single dielectric interfaces, we first show that the weak 
coupling (WC) Debye-Huckel (DH) theory remains quantitatively accurate up to the bulk ion density 
pt — 0.01 M, whereas the SC theory exhibits a good quantitative accuracy up to pt — 0.2 M, thus 
improving the accuracy of the DH theory by one order of magnitude in ionic strength. Furthermore, 
we compare the predictions of the SC theory with previous MC simulation data for charged dielectric 
interfaces and show that the proposed approaches can also accurately handle the correlation effects 
induced by the surface charge in a parameter regime where the mean-field (MF) result significantly 
deviates from the MC data. Then, we derive from the perturbative SC scheme the one-loop theory 
of asymmetrically partitioned salt systems around a dielectrically homogeneous charged surface. It 
is shown that correlation effects originate in these systems from a competition between the salt 
screening loss at the interface driving the ions to the bulk region, and the interfacial counterion 
screening excess attracting them towards t he surface. This competition can be quantified in terms 
of the characteristic surface charge a* — s/2p^JJjrlB), where Ib = 7 A is the Bjerrum length. In 
the case of weak surface charges a s <C at where counterions form a diffuse layer, the interfacial salt 
screening loss is the dominant effect. As a result, correlation effects decrease the MF density of both 
coions and counterions. With an increase of the surface charge towards <j*, the surface-attractive 
counterion screening excess starts to dominate, and correlation effects amplify in this regime the MF 
density of both type of ions. However, in the regime a s > <r* , the same counterion screening excess 
also results in a significant decrease of the electrostatic MF potential. This reduces in turn the 
MF counterion density far from the charged surface. We also show that for a s 3> o~t, electrostatic 
correlations result in a charge inversion effect. However, the electrostatic coupling regime where 
this phenomenon takes place should be verified with MC simulations since this parameter regime is 
located beyond the validity range of the one-loop theory. 

PACS numbers: 03.50.De,05.70.Np,87.16.D- 



I. INTRODUCTION 



The Poisson-Boltzmann (PB) formalism developed a 
century ago by Gouy [1] and Chapman [2] is still con- 
sidered today as the elemental theoretical description 
of electrostatic effects in various microscopic systems. 
The solution of the PB equations for charged macro- 
molecules immersed in salt solutions allows for example 
to determine protein folding pathways [3] or to under- 
stand the stability of DNA-binding proteins during their 
diffusion along DNA molecules One can also men- 
tion the microfluidic devices where the solute velocity is 
derived from the coupled solution of the PB and Navier- 
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Stokes equations [5]. Being a MF theory, the PB for- 
malism neglects however electrostatic correlation associ- 
ated with the interaction of the charged fluid with the 
system boundaries. It is thus clear that the PB equa- 
tion is a crude approximation for dielectrically discontin- 
uous systems such as water-air or water-membrane inter- 
faces where the electrolyte-surface interactions can signif- 
icantly exceed the thermal energy ksT in the proximity 
of the interface. 

The pioneering consideration of electrostatic corre- 
lation effects in inhomogeneous charged systems is 
doubtlessly Wagner's interpretation of the surface ten- 
sion excess of water with added salt in terms of the 
screened image charge interactions [6]. This theoretical 
framework that allowed Onsager and Samaras to derive 
their celebrated limiting law [7] was later improved in 
Ref. [8] by accounting for the non-uniform shielding of 
image interactions. Within a Wentzel-Kramers-Brillouin 
(WKB) approximation, the author ingeniously evaluated 
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the modification of the ionic self-energy and the surface 
tension by the interfacial variations of the ionic screening, 
which improved the agreement with experimental surface 
tension data. 

Correlation effects induced by the polarization charges 
at dielectric interfaces are also relevant to various indus- 
trial applications, among which one can mention wa- 
ter purification and desalination processes in artificial 
nanofiltration technology Non-linear electrostatic SC 
equations for confined electrolytes were derived in Rcf. [9 
within the Debye closure of the Bogoliubov-Born-Green- 
Kirkwood-Yvon (BBGKY) hierarchical equations. Ap- 
proximative solutions of these SC equations describe the 
selectivity of nanofiltration membranes in terms of a co- 
operation between the dielectric exclusion mechanism in- 
duced by image forces and the Donnan rejection driven 
by the membrane surface charge [10l [11] . They are fre- 
quently used today in order to predict experimental salt 
rejection rates in artificial nanofiltration processes |12j . 

The validity regime of the mathematical framework 
of these theories remained however unclear for several 
decades. The field theoretic formulation of the hetero- 
geneous Coulomb fluid model indeed provides some con- 
trol over the approximations involved in the consider- 
ation of correlation effects in nanoscale systems. The 
field has witnessed a dramatic growth during the last two 
decades. To give a non-exhaustive list, one can mention 
for example the consideration of electrostatic correlations 
in macromolecular forces for counterion liquids |13] and 
symmetric electrolytes [T3] at the gaussian level. The 
one-loop corrections to the density of counterions in con- 
tact with charged walls was also introduced in Ref . [TS] . 

Since the electrostatic coupling parameter increases 
with the ion density as T oc */p6, the exploration of the 
parameter regime beyond the dilute electrolyte limit ne- 
cessitates the consideration of non-linear effects neglected 
in gaussian field theories. As we will explicitly show in 
the present work, this is the regime where SC theories 
become relevant. The SC equations of Ref. [H] were red- 
erived in Ref. |16j within the field theoretic formulation 
of symmetric electrolytes. As stressed by the authors, 
these coupled SC equations are too complicated to be 
solved even numerically, and one has to make use of ad- 
ditional approximations in order to explore the underly- 
ing physics. These SC equations were solved in Ref. [T7] 
within a WKB-like approach in order to investigate ionic 
partitions around charged dielectric cylinders. However, 
we note that the approach used therein is not exactly 
a WKB method since the author did not make use of 
the WKB ansatz as in Ref. [8 , but rather solved the 
SC equations by assuming that the local screening pa- 
rameter in these equations does not vary with the spa- 
tial coordinate. A restricted variational theory was also 
proposed in Refs. [H] in order to understand non-linear 
effects in the process of dielectric exclusion from neutral 
slit nanopores. Furthermore, it was shown in Ref. |19j 
that the consideration of the interfacial ionic screening 
deficiency in this variational theory can considerably im- 



prove the agreement with MC simulations. An efficient 
restricted variational approach able to handle correlation 
effects induced by the surface charge from weak to strong 
coupling regime was also proposed in Ref. |2"D] . 

We introduced in Ref. |21j a simpler variational ap- 
proach for ions confined in charged slit nanopores that 
was able to handle the membrane charge with a good 
agreement with MC simulation results beyond the MF 
limit. We also applied this approach to cylindrical ion 
channels in order to show that the complications resulting 
from the curvature of the dielectric interface can explain 
the ionic current fluctuations in biological and artificial 
nanopores . We then extended the method by taking 
into account the excluded volume of ions in order to study 
excluded volume effects in the dielectric exclusion mech- 
anism [53] and macromolecular interactions [23]. Using 
similar ideas, we finally derived in Ref. |25j a non-mean- 
field dipolar PB equation in order to show that the in- 
terfacial solvent depletion at low dielectric surfaces can 
solely explain the low values of the experimental differ- 
ential capacitance data of carbon based materials. 

The weakness of the restricted variational approaches 
in Refs. [18-25] is that the restricted variational ansatz 
determines the nature of the final solution. The main 
goal of the present work is to overcome this limitation 
by solving the general SC equations of Ref. [TB] within 
two new computational approaches beyond the loop ex- 
pansion. The latter point will be shown to be crucial 
for understanding electrostatic correlation effects in di- 
electrically inhomogeneous systems where the one-loop 
theory is known to fail [55] . 

The article is organized as follows. We present in 
Section [II] an alternative derivation of the SC equations 
of Ref. [IB]. We then introduce in Section III two ap- 



proximative methods to solve these equations. The first 
method is based on an expansion of the formal inver- 
sion of the SC equations in powers of the fluctuating 
particle and charge density excesses around the weak- 
coupling theory. This approach is formally equivalent to 
the iterative solution of Hartree equations in condensed 
matter physics [57] ■ The second method is an exten- 
sion of the previously introduced WKB solution of these 
equations [5] at simple dielectric interfaces to the case 
of slit nanopores. We first compare in Section |IV| the 
predictions of these schemes for ion densities with the 
results of MC simulations that we ran for ions at neu- 
tral dielectric interfaces in order to establish the valid- 
ity domain of the WC and SC theories. Then, we test 
the validity of previous variational schemes for ions in 
slit nanopores, and also investigate within the WKB ap- 
proach the interaction between a charged rigid polymer 
and a dielectric wall. Furthermore, by comparisons with 
previous MC simulation data for ions in contact with a 
charged surface, we show that SC equations can handle 
correlation effects at charged dielectric interfaces beyond 
the WC parameter regime. Finally, we derive from the 
SC scheme the one-loop theory of asymmetrically dis- 
tributed salt solutions that we thoroughly examine. This 
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one-loop calculation bridges a gap between the DH the- 
ory of symmetric electrolytes at neutral interfaces |28j 
and the one-loop theory of counterfoil liquids in contact 
with charged interfaces [15] , Possible extensions of the 
concepts introduced in this article are discussed in the 
Conclusion part. 



II. REDERIVATION OF SC EQUATIONS FOR 
SYMMETRIC ELECTROLYTES 



We will derive in this part the self-consistent equations 
of Ref. [TH] in a shortcut way that does not require the 
evaluation of the Grand potential of the inhomogeneous 
electrolyte system. The grand canonical partition func- 
tion of a symmetric electrolyte composed of two species 
of valency ±q with q > 0, and fugacity Aj is given by 
a functional integral over a fluctuating electrostatic po- 
tential ^(r), Zq — JVcj) e~ H ^\ with the Hamiltonian 
functional [16] 



H[4>] = J dr 



[V<Kr)f 



— ia(r)(j>(r) 



2A,J dre E '- v - {r) cos [#(r) 



(1) 



where cr(r) is a fixed charge distribution, £b{ v ) — 
e 2 /(47rfcsTe(r)) the Bjerrum length, e(r) the static di- 
electric permittivity profile of the medium, and e the 
elementary charge. The wall potential V w (r) restricts 
the space volume accessible to ions. Furthermore, = 

2 

\v c {r — r')| r=r < is the self energy of ions in salt-free wa- 
ter, and v b c {r) = ^s/r the Coulomb potential in a bulk 
solvent, with £g = 7 A the Bjerrum length in a bulk elec- 
trolyte at ambient temperature T = 300 K. We finally 
note that in the present work, all energies are expressed 
in units of the thermal energy fc^T, the surface charge in 
units of the elementary charge e, and the dielectric per- 
mittivities in units of the dielectric permittivity of the 
air Eq. 

Our starting point is the compact form of the 
Schwinger-Dyson equation 



V(j) 



(r) 



-H[0]+/drJ(r)0(r) _ q 



(2) 



where J(r) is an external current. A rigorous proof of 
the equality ^ can be found for example in Ref. [25] . 
We will derive from Eq. |2]) two Ward identities relating 
the external electrostatic potential and the electrostatic 
propagator to higher order correlation functions. By act- 
ing now on the exponential with the functional derivative 
and setting J(r) = 0, one gets the following equation for 
the electrostatic potential, 



V£(r)V (<f>(r))+ia(r)-2Aiqe 



Ei-V w (r) 



sin [#(r)]> = 0. 
(3) 



We note that this equation was obtained in Ref. [301 . 
Furthermore, taking the functional derivative of Eq. (pi) 
with respect to J(r') and setting again J(r) = 0, we 
obtain a new relation 



k B T 



Ve(r)V<0(rMr')>+i<7(r) (0(r')) 



(4) 



-2A,ge 



Ei-V w {r) 



(0(r') sin [#(r) 



-5(r-r'). 



The approximation now consists in evaluating the aver- 
ages over the fluctuations in Eqs. Q and Q with the 
effective Hamiltonian of the most general quadratic de- 
pendence on the fluctuating potential 0(r) instead of the 
non- linear one in Eq. ([I]), 



i(f>o(r')/q] , 
(5) 



where the external electrostatic potential </>o( r ) 
— iq (4>(r)) and the inverse of the kernel i;(r, r') = 
(0(rWr')) - (0(r)) {<f)(r')) are solutions of Eqs. ^ 
and Q. The evaluation of the statistical averages in 
Eqs. Q and Q with the effective Hamiltonian Eq. ^ 
finally yields the self-consistent equations of Ref. [H] , 



Ve(r)V<fo(r) -e(rKe 



2 P -K,(r)- 



-Sv(r,r) ginh j^^] 



(G) 



e q 



a(v) 



Ve(r)Vu(r, r') - e(r) K 2 b e' v - {r) -^ Sv ^ r) cosh [0„(r)] v(i 
2 S(r-v% 



k R T 



(7) 



where we took into account the relation between the ion 

2 

fugacity and the bulk density A; = p bi e"^" Kt£f5 [T8 ] l2T ] . 
We also defined the ionic self-energy dressed with elec- 
trostatic correlations in the form 



Sv(r, r) = £ B Kb + v(r, r) - v^(0), 



(8) 



where we introduced the bulk screening parameter as 
kI = 8irq 2 £BPb- We also note that the local ion densities 
are given by [21] 



P±(r) = pue 



-V„(r)-i_^(r,r)T0o(r) 



(9) 



The relation Eq. ^ is a modified PB equation for the 
fluctuating external potential induced by the fixed sur- 
face charge around the MF potential. The second differ- 
ential equation ^ is a generalized Laplace equation that 
accounts for the local screening of the electrostatic propa- 
gator by mobile ions. In the next section, we will develop 
two approximative methods to solve these equations. 



III. COMPUTATIONAL SCHEMES 

We introduce in this section two computational 
schemes for solving the SC equations ^ and (JtJ) . The 
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first scheme consists in solving these equations by ex- 
panding their formal inversion around the one-loop elec- 
trostatic Green's function and the non-linear MF poten- 
tial in fluctuating excess charge and particle densities. 
The second scheme based on a WKB approach is an ex- 
tension of a previously introduced solution of the equa- 
tion Q for single neutral interfaces [H] to the more com- 
plicated case of neutral slit nanopores. The results of the 
theoretical approaches introduced in this section will be 
compared in the next section with extensive MC simu- 
lation data in order to determine the validity domain of 
Eqs. ((6) and Q. 



A. Perturbative solution of SC equations 

We present in this section an iterative solution of the 
SC equations ([6| and ^ around the MF external po- 
tential tp(r) and the one-loop Green's function vo(r,r'), 
respectively solutions of the equations [30 



where we introduced the expansion parameters A^, and 
A„ that will enable us to keep track of the perturbative 
order. These parameters will be set to unit at the end of 
the perturbative expansion. 

For charged liquids bounded by charged planar inter- 
faces located within the (x, y) plan, the translational 
symmetry considerably simplifies the problem. In this 
geometry, the external potential becomes simply a func- 
tion of the coordinate z, and the electrostatic Green's 
function can be expanded in 2D Fourier basis as 



v(r, r') 



4tt2 



e ik - r «i{z,z',k). 



(18) 



We now expand the Green's function and the fluctuating 
external potential in A„ and A,*, 



Ve(r)V<^(r) - e(r)^e~ v ' ro(r) sinh [<p(r)] 



2 

e q 



a(r) 
(10) 



v(r,r') = v (r,r')+ ]T X^5v nm (r,r') (19) 

n,m>0 

1>{z)= WVwOO- (20) 



n,m>0 



, Injecting these expansions into Eqs. ( 16 )-( 17 ) and 

|Ve(r)V - e(r)nle~ Vm( - r) cosh [<p(r)]j u (r, r') (11) Eqs. |l9])-(|2"o|), and keeping only the terms up to the or- 



5(r-r'). 



der A^A|, and \,A , we get 5v 01 (r 1 r / ) = 8v 02 (r,r') = 0, 
4>w(z) = if) 2 o(z) = 0, and 



The first step consists in injecting into Eq. ^ the fluctu- 
ating part of the external potential tp(r) = </>o( r ) — v( r )) 
and rearranging the resulting equation for ip( r ) with 
Eq. |7| in the form 

{ve(r)V - e{r)K 2 b e- v ™ (r) cosh [</>(r)]} ^i(r) (12) 
= e(r)n 2 b 6a{r) 

{ Ve(r)V - e(r)K 2 b e~ v ™ (r) cosh [^(r)] } v(r, r') (13) 

„2 

S(t - r') + e(r)nl5n(r)v(r, r'), 



k R T 



where we introduced respectively the fluctuating charge 
and particle density excesses as 



5a(r) = {n(r) sinh [</>o(r)] — sinh [ ( p( r )] 



(14) 



cosh [<p(r)] 4>(r)} 



,-V w {r) 



Sn(r) = {n(r) cosh [<fi (r)] — cosh [y(r)]} e 

(15) 

2 

with the Boltzmann factor n(r) = e~^ Sv ( r ' r \ Using now 
Eq. (11), the relations (12) and (13) can be inverted as 



5vxo(r,r') = -2p b q 2 J dziSn (zi)l2(r, r', zi) (21) 

Svn(r,r') = -2p b q 2 J dzin a ( Zl ) sinh [<p(zi]\ tpoi(zi) 
xI 2 (r,r', Zl ) (22) 

5v 2 o(r,r') = p b q A J dzin (zi) cosh [tp(z{)] 5v w (zi) 

xI 2 (r,r\ Zl ) (23) 
+ (2p b q 2 ) 2 J dzidz 2 Sno(zi)Sn (z2) 

xI 3 (r,r', Zl ,z 2 ) (24) 

for the corrections to the Green's function, and 

ipot(z) = ~2p b q 2 J dziv o (z,zi,0)Sa Q (zi) (25) 
i>n(z) = Pbq 4 / dziVo(z,Zi,0)no(zi)sinh[ip(zi)] 



x6v 10 (zi) 



(26) 



^(r) = -2p 6 g 2 A0 J driu (r, n)Sa(ri) (16) 
w(ry) = v (t,t') (17) 
-2p b q 2 X v J driUo(r,ri)^Ti(ri)u(ri,r'), 



^02(2) = ~2p b q 2 / dzi{j (z,zi,0)(5no(zi)V'oi(2 : i) 



(27) 



for the fluctuating potential. We note that we used above 
the notation Sv nm (z) = 5v nm (r, r) for the equal point ex- 
cess Green's function, and introduced the auxiliary func- 
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tions 

J 2 (r,r',zi) = 



47T 2 ' 



lk ' ( r " r !i ) v (z, z\ , k)v (zi , z' , k) 



(28) 



I 3 (t,t',zi,z 2 ) = 



in 2 ' 



where we have defined 
Pb = 



([ " lk e lk < r ^)vo(z,z 1 ,k)v (z 1 ,z 2 ,k) 



(35) 
(36) 



Xv (z 2 ,z',k), 



and 



Sa (z) = [n (z) - 1] sinh [<p{z)} e^™^ (30) 
5n (z) = [n (z) - 1] cosh [ip{z)\ e" y ™ (z) , (31) 

2 

with n (z) = e - 3 r 5vo ^\ Finally, the expansion of the 
ion densities in Eq. ^ yields at the same perturbative 
level the following expression. 



We also note that for ions confined in the nanoslit, the 
spatial integrals in Eqs. (21)-(23) run over the interval 
< z < d. Moreover, for a single neutral interface ob- 
tained from the slit pore geometry in the limit d — > oo, 
the weak-coupling Green's function ( 39 ) reduce to 



v (z,z',k) = 



2nt 



B I e ~Pb\z-z'\ 



Pb 



-Pb(z 



+ 2 ')}, (37) 



and the integrals in Eqs. (21 )-(23 1 have to be carried out 



P±{z) = pf\z)\l-K~5v w {z)T^m{z) (32) 



over the right half space z > occupied by the elec- 
trolyte. 



+X 2 v q 2 



? 2 

q q 

-y^n( z ) T ipu(z) ± — 5v w {z)tpoi(z) 



^Svj (z) - ^Sv 20 (z) 



7^01(2)^^02(2) 



where the zeroth order ion density is given by 

pf(z)=p bi e- v ^ z) -^ 5v ° { - z) ^ z \ 



2. Charged single interface 

We will derive in this part the zeroth order Green's 
function vv(z,z',k) and the associated ionic self en- 
ergy Svq(z) needed to compute the corrective terms in 
Eqs. (21 )-([2~7]) for the case of a dielectric interface lo- 
cated at z = 0, and carrying a negative surface charge 
a(z) = —<t s 5(z) with a s > 0. For this geometry, the MF 



potential solution of Eq. ( 10 1 is given by 



(33) 



We note that the above equations can be easily gener- 
alized to the case of asymmetrical electrolytes. We will 
derive below the one-loop propagator needed to compute 



tp(z) = 2 In 



1 



1 + e _Kf| ( z + 20 ) 



(38) 



the correction terms in Eqs. (21 \-b27h for the electrostatic 



Green's function and the external potential. 



1. Neutral interfaces 

In the case of neutral interfaces with a(z) — 0, the 
external potential and the corrective cross term of the 
Green's function both vanish, that is 4>q(z) = and 
Sv 11 = 0. Furthermore, we note that for electrolytes con- 
fined between two planar interfaces located at z = and 
z = d, and separating the solvent part with dielectric 
permittivity e w = 80 from the membrane matrix with 
permittivity e m < e w , the wall potential is defined as 
V w (z) = if < z < d, and V w (z) — 00 if z < or z > d. 

For the same slit geometry, the Fourier transform of the 
DH Green's function required to compute the correction 
terms in Eqs. (21|-(23l reads [3T] 



where we used the same notation as in Ref. [30] for 
the characteristic length of the interfacial counterion 
layer zq = — hx('j c (s)) / Kb , with the auxiliary parameter 
7c (s) = V s 2 + 1 — s, the dimensioneless parameter s = 
Ki,n, and the Gouy-Chapman length fi = l/(2irq£B<j s ). 

By injecting first the Fourier transform of the Green's 
function Eq. (18) into the Green's equation ( |l~T| ), we ob- 
tain 



—e{z)—v (z,z ,fc) 
az oz 



(39) 



-e(z)9(z) {p 2 b + 2«gcsch 2 [n b (z + z )}} v (z, z' , k) 

a 2 



k B T 



6{z-z'). 



(40) 



v (z, z\ k) 



2ir£ 
Pb 



For the ion source located in the water medium z' > 
0, the solution of this equation reads v(z, z' , k) — 
Cl e kz 6(-z) + [c 2 h+ (z) + c 3 h^ (z)] 9(z)9(z' - z), where the 
homogeneous solutions for z > were found in Ref. [30] 
in the form 



B 1 e ~Pb\z-z'\ 



(34) 



1 — A 2 .e~ 2pbd 



,-p b (z+z') , e p b (z+z'- 2 d) 



h ± (z)=e ±PbZ { IT — coth M 
I Pb 



Z + Zq) 



(41) 



The coefficients q for our system of impenetrable dielec- 
+2Abe~ 2pbd cosh (pb\z — z'Q] } , trie wall have to be computed by imposing the usual 
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boundary conditions associated with the continuity of the 
electrostatic potential and the displacement field [31) . 



v a (z = E_) =v(z = E+) 



s(z) 

dvp 
dz 



dv 



e(z) 



dv 



dz 



z— 

4tt£ b , 



(42) 
(43) 

(44) 



where E denotes the position of the interface at z = 0. 
One finds for < z < z' 

«(,(«, *',*:) = ^jf^ [M*)/l_CO + AM«)M2')] » 



(45) 



where the delta function is defined as 
A 



Kjjcsch 2 (K b z ) + (jp b - rjk) \p b - n b coth (K b z Q )} 

K^csch 2 (K b z ) + (p b + rjk) [p b + n b coth {n b z )] ' 

(46) 

with rj = e m /e w . In the limit a s — > where sur- 
face charge induced correlation effects vanish, the po- 
tential ( 45 ) is naturally reduced to the WC potential of 
Eq. (37). Finally, the electrostatic Green's function in 



real space follows from Eq. ( |18[ ) as 
dkk 



2tt 



Jo(A|r|,-r|,|)tJo(«y,*), (47) 



and the solution for z > z' is obtained by interchanging 
in Eq. (45) the variables z and z' . 



The computation of the corrections to the external po- 



tential in Eqs. (25)-(27) requires the evaluation of the 



function Vq(z, z', k) in the infrared limit, which reads for 
< z < z' 



v o (z,z',k^0) = — -e " 

K b 



{1 + coth [ Kb (z' + z )}} 
xH(z), (48) 
where we introduced the auxiliary function 



H(z) 



with 



K b z e 



A + ~il{s)n b z) e 



{< 



1 — K b z) e K 



A + 7c( s )( 1 + *bz) 



A = 



x coth [n b (z + z )) , 



2 7c 2 W 



T [7 c 4 ( S )-6 7c 2 ( S ) + l]. (51) 



Interchanging the variables z and z' in Eq. ( 48 ) yields 
the function Vo(z, z', k — > 0) for z > z' . 

Finally, the self-energy defined in Eq. ^ follows from 
Eq. ( 47 ) in the form 



[°° dk 

5v (z) = £ b kI / — r {-csch 2 [n b (z + z )] 

JO PbK 

Pb +co\h[K b {z + z Q )]\ e- 1 '""' 



(52) 



+A 



As expected, the potential Eq. ( 52 ) tends in the limit 
K b — > to the expression derived in Ref. [15] for the asym- 
metrically partitioned counterion-only system. Further- 
more, we note that the potential (52) is similar in form 



to the one derived in Ref. [30] for a symmetrical salt par- 
tition around a charged planar interface. However, due 
to the asymmetry of the salt distribution in our system 
as well as the presence of a dielectric discontinuity, the 
A function in Eq. ((52|) has a more complicated form and 



the potential (52) does not posses mirror symmetry with 



respect to the interface located at z — 0. Since we will 
exclusively investigate in this work ion densities in the 
water medium z > 0, we do not report here the expres- 
sion for the ionic self energy in the membrane medium 
z < 0. 



B. WKB approximation for neutral slit pores 



We will explain in this part the solution of the SC equa- 
tion ([7]) within a WKB approximation in the neutral pore 
limit <j{z) = where the external potential vanishes, i.e. 
4>{z) = 0. The calculation presented here is an exten- 
sion of a previous WKB solution of Eq. ^ for a single 
dielectric interface [8] to the more complicated case of a 
slit nanopore. Unlike the approach of Section [ill A| based 
on a perturbative expansion of the formal solution of SC 
equations, the WKB approximation has the advantage of 
being a non-perturbative approach. Consequently, while 
moving from the bulk towards the pore wall, it will be 
shown below that the WKB solution can interpolate be- 
tween the DH regime with k{z — > 00) = n b and the dilute 
electrolyte regime with n(z — > 0) = in a self-consistent 
way. 



First of all, by injecting the Fourier expansion Eq. ( 18 ) 
into Eq. Q, one gets 



(49) 



[~d z e(z)d z + e(z)p 2 (z)] v(z, z , k) 



k B T 



S(z- 



(53) 



(50) where we introduced the function p{z) = ^/k 2 + k 2 (z) 
whose fc-dependence is implicit. Furthermore, we defined 
a local screening parameter in the form 



K 2 (Z) 



K h e 



-V m (z)-S T Sv(T,r) 



(54) 



In the present work, we will exclusively need the Green's 
function associated with source particles within the slab, 
i.e. < z' < d. Without making any approximation yet, 
the general solution to the second order differential equa- 
tion ( 53 ) in the slit geometry can be expressed in terms 



of its two independent homogeneous solutions h±{z) in 
the form 
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v(z, z', k) = Cl e kz 0(-z) + [c 2 h + (z) + c 3 h_(z)] 9{z)9{z' '- z) + [c 4 h+{z) + c 5 h_(z)} e(d-z)e{z-z , )+c 6 e- kz e(z-d), (55) 



where the coefficients Ci are integration constants to be found by imposing the boundary conditions |42j)-(44) at the 
interfaces £ = and £ = d. After some long algebra, we get 



v(r,r')=£ B 
where 



dkk 



G(z') 
F(z') 



J-h (fc|r,| -r||]) {[e w h'_(0) - e m kh-{0)] h+{z) - [e w h' + (0) - e m kh+(0)] h^(z)} , (56) 



G(z)=e w [h+(z)h'_(d) - h_(z)ti + (d)] + e m k[h + {z)h-{d) - h_{z)h + {d)] 



(57) 



and 



F(z) = [h-(z)h' + (z) - h + {z)ti_{z)} {e 2 w [h' + {Q)h'_{d) - ti_{ti)ti + {d)} + e 2 m k 2 [h-(0)h + (d) - h + (0)h-(d)] (58) 

+e w e m k [hf + (0)h-(d) - ti_{d)h + {0) + h_{Q)h' + {d) - h + (d)ti_(0)] } . 



The formal solution ( 56 ) can be actually used in order to 
solve Eq. |7|) exactly by evaluating the functions h±(z) 
from a numerical solution of Eqs. (53) and (54) in a two 



dimensional lattice (k, z) by iteration. We will however 
leave the exploration of this idea for a future work, and 



evaluate Eq. ( 56 ) within the WKB approximation where 



the homogeneous solutions of Eq. ( 53 1 read [5] 



h±{z) = [p{z)] 



-1/2 



exp 



± / dz'p(z') 
o 



(59) 



Furthermore, in order to simplify the numerical task, we 
will limit ourselves to the case e m < e w characterized by 
a continuously varying k(z) on the interval < z < d and 
vanishing ion densities on both sides of the interfaces, i.e. 



lim p(z) = lim p(z) = 0. 

2^0± z^d± 



(60) 



After some long but straightforward algebra, one obtains 
the kernel in the form 



v(r, r') = i; 6 (r, r') + v d (r, r'), 
where the bulk part is given by 

dkk 



(61) 



v b (r,r') = it 



= J 



'o ^/p{z)p(z') 

Xe -\H*,*')\ 

and the dielectric contribution reads 

dkk J o( fc l r ll - r '|l) 



(62) 



Ud(r,r') = ^sAo 



x ( e -I(0,z)-I(0,z') + e -I(z,d)-I(z\ 

+2A e- 2/(0 ' d) cosh [I(z,z')]\ 



(63) 



Wc defined in Eqs. (62) and (63) the integral 



I(zi, zi) = / dzp(z) 

Jz! 

and introduced the dielectric jump function 



Ao 



(64) 



(65) 



It is interesting to note that the bulk part of the elec- 
trostatic Green's function Eq. ( 62 ) is not of the DH form 
i>£>#(r,r') = £b lre~ KbT ', with r = |r — r'| the interi- 
onic distance. The difference between both potentials is 
clearly due to the breaking of the spherical symmetry by 
the non-uniform ionic screening, an intcrfacial effect ab- 
sent in the DH theory. Moreover, Eq. (62) shows that by 



approaching the dielectric surface where the ion density 
vanishes, the potential (62) tends to the usual Coulomb 
law, i.e. lim z 2 /_>.o+ Vb(r, r') = £b/t- In the next part, we 
will take advantage of the ability of the WKB solution to 
self-consistently interpolate between the dilute limit and 
the DH regime in order to evaluate the asymptotic small 
distance limit of polymer-interface interactions. 

For the computation of the local ion densities, we ex- 
clusively need the self energy of ions defined in Eq. (JsJ , 
which reads 



Sv(r, r) = Sv(z) = Ib [nb — k(z)] (66) 
dkk e - 2/ (°' z ) + e - 2I{ - z ^ + 2A n e- 2/ (°- d ) 
p(z) 



1 - A 2 e-2/(o,<i) 



We notice that the first term on the rhs of this equation 
is associated with the local variation of the ionic cloud 
around the source ion, and since n(z) < kj, for e m < e w , 
this solvation energy positively adds to the image charge 
repulsion contribution, i.e. the second term on the rhs of 
Eq. ([66]). 

The relations ( 54 ) and ( [66] ) form a set of closure equa- 
tions that has to be solved by iteration on a bidimensional 
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FIG. 1: (Color online) Ion density profiles at the dielectric interface against the distance from the surface with e m = 1, e w = 80, 
and ion diameter a,i — 4.25 A at the bulk ion concentration (a) pb = 0.01 M, (b) pb = 0.1 M, (c) pb = 0.2 M, and (d) pb = 0.4 M. 
The red lines are our MC simulation data, the blue lines the WC theory, and the dashed brown and black lines are respectively 
the WKB and the second order perturbative solutions of the SC equation |7|. The black dots in (a) denote the iterative solution 
of the closure relations Eqs. (541 and (66 1 explained in Appendix [A} while the blue squares in (b) are the simulation data from 
Ref. [33]. The dashed blue line in (c) marks the third order perturbative solution of SC equations. 



lattice (k,z). The iterative approach consists in inject- 
ing into Eq. ( |54[ ) the ionic self energy obtained from the 
WC potential Eq (39), which yields an updated screening 
function K\{z). One then evaluates the new self energy 



profile from Eq. (66) with k,\{z), and the iterative cy- 



cle continues until self-consistency is achieved. We note 
that this iteration method was used in Ref. [5] to com- 
pute an analytical expression for the ionic self energy at 
single dielectric interfaces with e m — 0. An extension 
of this calculation to finite values of e m is presented in 
Appendix [X] Before concluding, we also note that for 
simple interfaces, the ionic self energy ( 66 ) reduces to 
the expression derived in Ref. [5] , 



5v{z) =l B [n b - k{z)]+£ b A 



dkk 

PcO) 



-2/(0, z) 



(67) 



IV. RESULTS 



Neutral interfaces 



We will first determine in this part the salt concen- 
tration range where the SC equation (|7|) remains quanti- 
tatively accurate for electrolyte systems in contact with 
neutral dielectric interfaces separating the solvent part 
with ions from a membrane region free of ions. This ge- 
ometry is relevant to the water-air interface as well as to 
membrane nanopores characterized by strong ionic con- 
finement effects. The computation schemes developed in 
the previous parts will be then applied within this range 
to slit nanopores and polymer-interface systems in order 
to test the validity of the DH theory and a restricted vari- 
ational approach [21] relevant to experimental nanofiltra- 
tion studies 1121. 
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this conclusion at an analytical level. Indeed, the first 
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FIG. 2: (Color online) Pore averaged ion densities in slit pores 
against the pore size for pt = 0.1 M, e m = 1, and e m = 80. 
The light blue line is the WC theory, the black line is the 
SC theory, the red line is the WKB solution, and the dark 
blue line denotes the restricted variational scheme of Ref. [21] . 
The inset displays local ion densities in the pore for the same 
model parameters. 



1. Ion densities 

We establish in this part the validity domain of the 
SC equation (J7j) for neutral single interfaces and slit 
nanopores by comparing the predictions of the theoretical 
schemes developed in the previous parts with MC simu- 
lation data for ion densities. The details of our numerical 
simulations can be found in Appendix[Cj We note that in 
order to be able to compare our simulation results with 
previous MC simulation data from Ref. [33J , we chose the 
diameter of ions as a.; = 4.25 A, and the dielectric per- 
mittivity of the membrane and the water respectively as 
s m = 1 and e w = 80. The comparison of our MC data 
in Fig. [ijb) with the numerical results of Ref. [33] shows 
the good agreement between the two simulation results. 

The comparison of the WC theory in Fig.[TJa) with MC 
data shows that the concentration pb = 0.01 M marks 
the boundary of the ion density range where the WC the- 
ory starts to deviate from the simulation result, although 
even at this concentration, it exhibits a reasonably good 
quantitative agreement with the MC data and the SC 
theory. We note that this density corresponds to an elec- 
trostatic coupling parameter T — k^b — 0.2. For larger 
bulk concentrations in Figs. [TJb-d where this deviation 
becomes more pronounced, one sees that the WC theory 
systematically overestimates the ion density. Within the 
restricted variational theory of Refs. [TBI EH]> this over- 
estimation was shown to originate from the unability of 
the WC image potential of Eq. (37) to account for the 
reduction of the ionic screening at the interface. The 
transparency of the present method allows to confirm 



shows with Eq. (31) that due to the interfacial ion defi- 



ciency Suq(z) < 0, the positive correction to the image 
potential 6vio(z) > increases the amplitude of the WC 
potential, leading to a stronger dielectric exclusion at the 
interface. 

An inspection of Figs. [TJb-d shows that the predic- 
tions of the WKB and the perturbative solutions of the 
SC equation ^ both exhibit a good agreement with the 
simulation data up to pb = 0.2 M, thus improving the 
quantitative accuracy of the WC theory approximately 
by one order of magnitude in ionic strength. Hence, the 
deviation of the SC theory from the simulation data tak- 
ing place at pb = 0.2 M (or T ~ 1.0) establishes this value 
as the characteristic density where the quantitative ac- 
curacy of the SC theory breaks down. We note however 
that it is unclear whether the failure results from electro- 
static correlation effects, or excluded volume effects not 
included in the SC theory that start to set on. This point 
can be enlightened in future by solving the extended SC 
equations of Ref. |23j that can account for ionic excluded 
volume effect. 

Before concluding the discussion of the ionic partition 
at single dielectric interfaces, we would like to note two 



points. First of all, we show in Fig. [I 



that the approxima- 



tive solution of the closure relations ( 54 ) and ( 66|) by iter- 
ation explained in Appendix|X]fits very well tl 
ical solution in the dilute regime pb 



;heir numer- 
0.01 M. This ob- 
servation is useful since this approximative solution will 
be used in the next section for an analytical evaluation 
of the interaction of a rigid polymer with the dielectric 
interface. Then, we emphasize that the perturbative SC 
solutions reported in Fig. [TJ are obtained at the second 
order perturbative level introduced in Sec. Ill A To as- 
certain the convergence of our perturbative scheme, we 
also reported in Fig. JlVc) the third order calculation ex- 
plained in Appendix IB] (see the dashed dark blue curve). 
A careful inspection shows that the third order result 
is hardly distinguishable from the second order result, 
which confirms that for neutral single dielectric inter- 
faces, the second order perturbative solution is sufficient 
within the validity domain of the SC equation ([7]). 

We now illustrate in Fig. [2] the ion densities and par- 
tition functions k = (p{z)/ ' Pb) p in slit pores with e m — 1 
and e w = 80, in contact with a bulk reservoir of ionic 
concentration pb = 0.1 M. The plots compare the re- 
sults obtained from the WC theory, the perturbative 
and WKB solutions of the SC approach, and the gener- 
alized Onsager-Samaras (GOS) approach introduced in 
Ref. [3TJ. We note that the model parameters in this 
figure were chosen in such a way that the results ob- 
tained from the second and the third order perturbative 
solutions of the SC equation ^ are practically superim- 
posed, thus guaranteeing the convergence of the pertur- 
bative solution. An inspection of the main plot and the 
inset of Fig. [2] shows that as in the single interface case, 
the WC theory underestimates the dielectric exclusion, 
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while the WKB and the perturbative solutions of the SC 
theory are again very close to each other. Moreover, one 
notices in the inset of the same figure that the GOS for- 
malism is able to accurately reproduce the ion density 
close to the interface, but slightly overestimates the di- 
electric exclusion in the mid-pore area. However, we see 
in the main plot that despite this weak discrepancy, the 
GOS approach remains very accurate in estimating the 
partition coefficients over a large interval of pore thick- 
ness. This observation might indeed explain the success 
of GOS-like solutions of the SC equation ^ frequently 
used in artificial nanofiltration studies in order to esti- 
mate exprimental salt rejection rates [TUHI2- 

2. Polymer-surface interactions 




O05 (15 5 



In this section, we will evaluate within the WKB ap- 



proach introduced in Section HI B the energetic cost to 
drive a rigid polymer from the bulk reservoir to the prox- 
imity of the interface. If the surface charge density of the 
polymer r is weak enough so that it does not significantly 
affect the interfacial ion densities, the electrostatic energy 
of the polymer located at (y = 0, z) can be obtained from 
the relation [2H1 



F(z) 



dridr 2 



<r(ri)?;(ri,r2)c7(r2), 



(68) 



where the electrostatic Green's function is given by 
Eq. (61 1, and the linear charge density of the polymer 
is 



a(r') = TS(y')S(z' 



(69) 



The net energetic cost for bringing the polymer from the 
bulk to the dielectric surface located at z = is given by 



A/(z) = i [F(z) - F b ] 



(70) 



where L is the polymer length, and the total electro- 
static energy of the polymer in the bulk electrolyte reads 
F b = lim^oo limbec F(z). We note that within the DH 
theory, the same energy density was derived in Ref. |28j 
in the form 



Ibt 2 



dk y 

Pb 



A b e- 2p » z , 



(71) 



where p b and Af, are respectively given by Eqs. (35) 



and (36 1. In the asymptotic limit of small polymer sur- 



face separations, Eq. (|71j) was also shown to reduce to a 

(72) 



simple logarithmic law, 

Af DH (z) ~ A„e B T 2 ln(l/K b z) 



Evaluating the integrals in Eq. (68) with the Green's 
function Eq. (61) in the limit d — ¥ oo and the polymer 



FIG. 3: (Color online) Reduced polymer-interface free energy 
against the distance from the dielectric surface for the mem- 
brane permittivity e m = 1, and salt concentrations pt, = 0.01 
(red curves) and 0.1 M (blue curves). The solid lines are from 
the WKB theory and the dashed lines denote the prediction 
of the DH theory. The black dashed lines mark for each curve 
the asymptotic small z behavior. 



charge distribution of Eq. ( 69 ) , we obtain the free energy 
profile as 



A/(z) =£ b t 2 In [ Kb / K (z)}+£ B T : 



■A 



dfe 

-oo P{z) 



V e -2/(0,z)^ 



(73) 

We emphasize that since the WKB solution was derived 
in Sec. |IIIB"| for the case of a dielectric discontinuity where 
the ion density vanishes at the interface, the relation ( 73 ) 
is valid in the permittivity range e m < e w . 

It is seen that Eq. (73) is composed of a local ionic 



solvation part absent at the DH level, and a second term 
associated with the dielectric jump at the interface. The 
solvation term that depends logarithmically on the ionic 
density accounts for the spatial variations of the ionic 
cloud density around the charged polymer strand. Since 
the screening of the polymer charge lowers its free energy, 
this contribution acts as an additional force pushing the 
polymer towards the bulk area where the ionic density is 
maximum. 

DH and SC free energy profiles of Eq. ( 71 1 and ( 73 ) 
are illustrated in Fig. [3] for the membrane permittivity 
e m = 1, and salt concentrations p b — 0.01 and 0.1 M. We 
first notice that the free energy barrier evaluated within 
the SC theory is significantly larger than the prediction 
of the WC theory. Then, one notices that in the close 
neighborhood of the interface, the same barrier is char- 
acterized by a regime that is clearly independent of the 
bulk salt concentration. 

In order to elucidate these two points, we will evaluate 
the asymptotic small z limit of Eq. ( 73 ) within the pertur- 
bative approach explained in Appendix [X] To this end, 
we inject into Eq. (73) the screening function Eq. (54) 
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evaluated at the first order iterative level with the den- 



sity profile Eq. ( A5 ) , 



where we defined 

Kq(z) = K^exp 



l-\5v x {z) 



q 2 t B A 2Kb 

Az 



(74) 



(75) 



and expand the result in powers of Svi(z). At the leading 
order O (Svi(z)), we get 



A/(z) ~ e B r 2 ln[ Kb /K (z)}+e B T 2 ^ 



dk 



y_ e -2i (o,z) 



5vi(z) 



(76) 



A Q £ B q 2 T 2 



dk 



o Po(z) {2Po{z) 

,K 2 (Z>) 



+ / dz'^- L 5v 1 (z')\e 
J a Po{z') J 



-2/ (0,a) 



where the subscript means that the screening function 
k(z) in the functions p(z) and 1(0, z) should be replaced 
with Kq(z). We will now derive from Eq. ( 76 1 the asymp- 
totic limit of the free energy close to the interface. If we 
note that the screening function n[z) is very small with 
respect to K b in the neighborhood of the surface, we can 
assume that Po(z) — \/k 2 + n 2 (z) varies slowly close to 
the dielectric interface. One can thus take the function 
Po{z) out of the integral in Iq(0, z), and carry out the in- 
tegration over k y . Expanding the result up to the linear 
order in z. one obtains 



A/(z) 



a (i + a k!tV 



-"/A £ b t 



+ A q £ b t 2 ln(l/K b z) 



U%r\ b [ 1 



rA? 



(77) 



where 7 ~ 0.5772 stands for the Euler's constant. The 
asymptotic law ( 77 1 reported in Fig. [3] by dashed black 
lines is shown to accurately reproduce the behavior of 
the numerically computed free energy barriers close to 
the surface. We see in Eq. (77) as well as in Fig. [3] that 
the logarithmic dependence of the energy barrier on z 
predicted by the WC theory is actually dominated by 
an unscreened algebraic decay. This algebraic regime as- 
sociated with the interfacial dielectric exclusion of ions 
explains the salt free portion of the free energy barrier as 
well as its strong amplitude observed in Fig. [3j 



B. Charged interfaces 

We apply in this part the perturbative SC scheme in- 
troduced in section |III A| to charged single interface sys- 
tems in contact with a symmetric electrolyte composed 




■ MC 

Mean-field 

Self-consistent 



4 , 6 

z/a. 



(a) 



(b) 




Cl 




(c) 

FIG. 4: (Color online) (a) Electrostatic potential profile for 
Pb = 0.01 M, and counterion (upper black curves) and coion 
(lower red curves) densities for pt = 0.01 M (b) and 0.1 M (c). 
The solid curves are the SC theory, the dashed lines denote 
the MF result, and the squares are the MC simulation data 
from Ref. [34] with ion diameter ai = 4.25 A, e m ~ 1, and 
s w = 78.5. 



of monovalent ions q = 1. By comparing the theoreti- 
cal predictions for ion densities and the interfacial elec- 
trostatic potential profile with MC simulation data, we 
will first investigate electrostatic correlation effects at 
charged dielectric interfaces. The perturbative SC ap- 
proach presents itself as a useful computational tool par- 
ticularly in this case of a dielectrically discontinuous sys- 
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tern where the one-loop expansion is known to fail |26j . 
We will then expand in the limit s m = e w the perturba- 
tive solutions of the SC equations ^ and Q introduced 
in section |III A| at one-loop order, which will provide us 
with an analytical one-loop theory of asymmetrically par- 
titioned symmetric electrolyte solutions around charged 
planar interfaces. 



Correlation effects at charged dielectric interfaces 



We illustrate in Fig. |4^a) and (b) the prediction of the 
MF theory and the second order perturbative solution 
of the SC theory for the electrostatic potential and the 
ion densities close to a charged interface with the surface 
charge a s — 5.5 x 10~ 2 e nm~ 2 , the bulk ion density pb — 
0.01 M, and the dielectric permittivity values e m = 1 and 
e w = 78.5. Furthermore, Fig. |i|c) displays the local ion 
densities for the same dielectric constants, but with a 
slightly stronger surface charge a s = 7.75 x 10 -2 e nm -2 
and a significantly higher salt concentration pb = 0.1 M. 
We also reported for the three plots the MC simulation 
data from Ref. [34] for ions of diameter a% — 4.25 A. 

The first point to be noted in these plots is the good 
quantitative accuracy of the SC theory in predicting the 
electrostatic potential profile and the ion densities ob- 
tained from MC simulations in the neighborhood of the 
charged dielectric interface. This is particularly notice- 
able for the case pb = 0.1 M where the MF result ex- 
hibits a clear disagreement with the simulation data. 
The deviation of the MF curves from the MC data can 
be explained as follows. First of all, the equations (25) 



and ( 30 1 show that the ionic screening deficiency induced 



by the dielectric exclusion at the interface gives rise to 
a positive charge density excess, i.e. 5<j s {z) < 0, which 
in turn results in a negative correction ij}Qi(z) < to the 
MF potential, as can be seen in Fig. |4]ja). Figs. |4|b) 
and Qc) show that the underestimation of the strength 
of the electrostatic potential by the MF theory is respon- 
sible for an underestimation of the counterion attraction 
and the coion repulsion by the surface charge for z > 2 
A. This means that in the presence of a strong dielectric 
discontinuity between the solvent and the weakly charged 
substrate, correlation effects reduce the amplitude of the 
MF level charge separation. Moreover, because the MF 
theory is unable to account for the charge-image repul- 
sion that dominates the ion-surface charge interaction for 
z < 2 A, both coion and counterion densities are over- 
estimated by the MF theory in the close neighborhood 
of the interface. However, it will be shown in the next 
part that in the case e m = e w , this picture is gradually 
reversed with increasing surface charge. 



2. One-loop theory of asymmetrically partitioned 
electrolytes 

We present in this part the one-loop theory of asym- 
metrically distributed electrolytes around a charged sur- 
face that will be shown to directly follow from the pertur- 
bative SC scheme developed in section |III A| This one- 
loop calculation that will allow an analytical investiga- 
tion of electrostatic correlation effects bridges a gap be- 
tween the DH theory of ions at neutral dielectric inter- 
faces [28] and the one-loop theory of counterion liquids 
at charged interfaces [15] . 

By rescaling first all lengths with the screening pa- 
rameter according to z = KbZ in Eqs. (21)-(27) and in 
Eq. (52), and expanding in the limit e m — e w the same 



equations in powers of the electrostatic coupling param- 
eter r = iB^b, one obtains the following series 

Sv (z) = TSv {z;s) (78) 

Sv nm (z) = rS ~ v nl&s) (79) 



z>n+m+l 



1pnm(z) = E ^1(1; S ), 



(80) 



z>n+m 



where the functions under the bar sign depend exclusively 
on the dimensionless distance z and the parameter s = 
KbH- We note that for e m ^ e w , the self-energy 8vq(z) 



defined in Eq. (52) and appearing in the argument of the 
exponentials in Eqs. (30 l-plj) becomes singular at z = 0. 
Thus, such a loop expansion is valid exclusively in the 
absence of a dielectric discontinuity [26] . 



The equations (78)-(80) show that the only potentials 



surviving at the one-loop level are the ion self-energy 
Svq(z), which is purely on the order T, and the one-loop 
part of the correction V'oi(z) to the MF external potential 
ipu(z) = r-0Q^(z; s), which in turn reads 



ipu{z) = p b q / dz 1 v (z,zi,0)8vo(zi)smh[ip(z 1 )]. 
Jo 

( 81 ) 

A simpler expression for this one-loop potential correc- 
tion will be given below. Furthermore, by expanding 
Eq. (32) up to the order 0(T), the one-loop density fol- 



lows in the form 



p±{z) = pue 



1 - y*«o(2) T ipil(z) 



(82) 



Expanding now the differential equation ( 12 ) satisfied by 



the fluctuating part of the external potential ip(z) at the 
one-loop order, one obtains 



d 2 iPu(z) 



dz 1 



cosh [ v {z)]i> u {z) 



- v ^Sv {z) sinh [<p(z)] 



(83) 



The one-loop correction to the external potential was 
computed in Ref. [30] by solving the equivalent of the 
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equation ( |83| ) for the charged interface system of a sym- 
metric salt distribution, which is formally equivalent to 
the calculation above since we note that the inversion 
of Eq. (83) directly yields the expression (81). Fur- 
thermore, by integrating Eq. (83) from z = — oo to 
z = +oo and using the expression (81), one can show 
as in Ref. [30 that the one-loop solution (81) automati- 
cally satisfies the global electroneutrality condition, that 
is qf Q dz[pi+(z)— P-(z)\ — (T s . In other words, as 
stressed in Ref. [15] for the case of the inhomogeneous 
counterion liquid in contact with a charged interface, the 
total charge density gets a vanishing contribution from 
the one-loop theory. Moreover, we note that the one-loop 
potential ip{z) + ipn(z) satisfies the Gauss law since one 
finds from Eq. (81) that the one-loop correction has a 
vanishing derivative on the surface, i.e. V'i;(0 + ) = 0. 

We now note that performing the variable transforma- 
tion k — y u = Pb/nb in the integral of Eq. (52), the ex- 
pression for the ionic self-energy can be recast in a more 
manageable form 



5vq(z) = r 



^-{-csch 2 [z-ln 7c ( S )] (84) 

XL 1 — 1 1 



+A (u + coth [z — ln7 c (s)]) e 
where we introduced the function 



1 + (u - Vu 2 - 1) s (su - Vs 2 + 1) 

A — ; , : ; . - ■ (85) 

1 + (u + Vu 2 - 1) S (su + y/s^Tl) 

In the regime of weak surface charges or high salt concen- 
trations corresponding to a large value of s, the integral 
in Eq. (84) can be analytically evaluated by expanding 



the integrand in powers of 1/s. At the leading order, one 
finds 



Sv (z) = Sv 0n (z) + s 2 Svq c (z) + 0(s 4 ) 



(86) 



where the vanishing surface charge part was previously 
derived in Ref. |21j in the form 



Sv 0n (z) (I + Z) 2 _r 

— -e 



1 



K 2 (2z) 



(87) 



and the surface charge contribution is now given by 

Svo c (z) 



1 



1 



-25 



K (2z) 



1 , 



l + z(z--)(l + e- 22 ) 



Ki(2z) 
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3 




[h 


- 4 
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" 2I 2 


+- 
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-42 _ 




2z A 







ln(4z) 



Ei(-4z). 



-1z 



The above equations make use of the modified Bessel 
functions of the second kind K n (x) and the exponential 
integral function Ei(a;) |35j . The close form expression 
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FIG. 5: (Color online) (a) Ionic self energy renormalized by 
the coupling parameter Y against the reduced distance from 
the interface KbZ for different values of the parameter s. See 
text for details. 



2 with the integral form of 
by 



Eq. ( 86 1 is compared for s 
Eq. Q in Fig. [5j 

First of all, we note that the potential given 
Eq. (87) originates from the ion deficiency in the mem- 
brane medium z < 0, and it is known to bring a purely 
repulsive contribution to the potential Svq(z) |21j . This 
repulsive part of the potential marks the upper bound- 
ary of the self energy curves in Fig. [5j The second part 
5vo c (z) of Eq. ( 88 ) induced by the surface charge is purely 
negative. As shown in Fig. [5j this term brings in turn 
a net attractive contribution to the self energy. In the 
asymmetrically distributed salt system considered in the 
present work, it will be shown that electrostatic correla- 
tion effects are mainly driven by the competition between 
these two opposite mechanisms. This competition will be 
thoroughly investigated below. 

For small separations from the surface z«l, Eq. (88) 
takes the asymptotic form 



Svp(z) 

r 



[4 7 -3 + 41n(z)] 



(89) 



1 



3s 2 
+0(z 2 ) 



|l2 ln(2) - 7 - | [5 + 127 + 12 H^)} } 



The equation (89) shows that the potential 5vq(z) ex- 
hibits a linear decay with the distance z close to the in- 
terface. We recognize in the first term of this expression 
the energetic cost £bK{,/3 to drive an ion from the bulk 
to a neutral interface separating a membrane medium 
free of ions and a salt solution of ionic strength k 2 [3T] . 
Moreover, the contribution from the surface charge is 
shown to lower this barrier at the interface by the amount 
-(121n2-7)r/(3s 2 ) ~ -0.44^ s /(k^ 2 ), and also to give 
rise to a potential minimum located between the maxi- 
mum counterion concentration at z = and the bulk 
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region z — oo. By comparing the magnitude of these two 
terms, one finds that at the particular value /i -1 ~ Kb 
where the thickness of the interfacial counterion layer p 
becomes equal to the radius of the ionic cloud 1 around 
a central charge in the bulk area, the potential vanishes 
on the surface, 6vo(0) = 0. Furthermore, the self energy 
profile 5vq(z) remains purely attractive for higher sur- 
face charge values or lower bulk ion concentrations. This 
point is also illustrated in Fig. [5| 

We note that one can also express the competition be- 
tween the bulk and interfacial solvation effects in terms 
of the density of the interfacial counterion layer c ~ 
I-kIboI [IS] and the bulk salt density as 4pb < c. This 
inequality characterizes the regime where the solvation of 
ions by the interfacial counterion layer attracting them 
towards the surface overcomes the strength of the interfa- 
cial salt screening loss driving the charges far away from 
the interface. We finally note that the parameter do- 
main where the potential 5vq(z) is negative can be also 
rewritten in terms of the balance between the bulk ion 
concentration and the surface charge as pb < 7t^bc 2 /2. 
It is interesting to note that this inequality is indepen- 
dent of the ion valency. For instance, for a dilute salt 
solution with concentration pb = 0.01 M, this inequality 
becomes an equality for the characteristic surface charge 
ct s = 7.5 x 10~ 2 e nm~ 2 . 

In the opposite regime of large separations from the 
surface z» 1, the potential Eq. (86) takes the asymp- 
totic form 



Sv (z) 



-2z 



2z 



1 



-4z\ 



+ 7Z + z ln(4z) 



(90) 



In the strict limit z — > oo, the leading term of this asymp- 
totic law reads 8vq(z)/T ~ — s -2 ln(4z)e~ 2z . This shows 
that far enough from the interface, the attractive surface 
charge contribution will always dominate the repulsive 
solvation force associated with the salt screening loss in 
the proximity of the interface. 

In the opposite small s (or the Gouy-Chapman) 
regime, Eq. ([84l yields 
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FIG. 6: (Color online) (a) One-loop correction to the external 
potential Eq. ( |95[ ) renormalized by the coupling parameter T 
against the dimensioneless distance KbZ. The dashed black 
lines denote the large distance asymptotic limit Eq. ( 99 1. The 



inset displays the surface potential and the renormalization 
factor behind Eq. (99 1 against the parameter s. 



The large distance asymptotic limit z ^> 1 of Eq. ( 84 ) 



can be computed for an arbitrary finite value of s in the 
form 



Sv (z) 



2^(s)e- 2z [e 42 Ei(-4z) - ln(42) - 7 J (94) 
+0 (e~ 42 ) . 



This asymptotic law is reported in Fig. [5] for 8=1. Not- 
ing that lim^oo e x Ei(—x) — —1/x, the strict large dis- 
tance limit z — > oo of Eq. ( |94| is obtained as 5vq(z) ~ 
— 2Tj^(s)e~ 2z ln(4z), which is a purely negative function. 
One can verify that this equality consistently recovers the 



Sv (z) 



1 

2I 



leading terms of Eq. (90 1 and Eq. (93) in the limits of 



{e 2z — [7 + ln(4z) — Ei(— 4z)] z csch (2)} large and small s, respectively. The relation (94) con 



-0(8). 



(91) 



In Fig. [5J it is shown that the limiting law Eq. ( |9"lj ) is 
purely negative. One also notices that this asymptotic 
limit marks the lower boundary of the self energy curves. 
For small separations from the surface z <§C 1, the poten- 
tial decays algebraically with increasing distance, 



firms the conclusion that we previously reached for weak 
surface charges : in the presence of an arbitrary finite sur- 
face charge and salt concentration, the ionic self-energy 
will always posses an attractive branch far enough from 
the interface. 

We will now investigate the behaviour of the one-loop 
correction to the external potential ipu(z) with respect to 
the parameter s. Evaluating first the integral in Eq. (81 ) 



( 92 ) with Eqs. d3Sb , ((481), and ([52]), one gets for the one- loop 



and for large separations z 3> 1, it is screened exponen- 
tially, 

~ -2 [7 + ln(4z)] e- 2z ~ + O (e^) . (93) 



correction to the external potential 



i>il(z) = — rcsch [z - ln 7c (s)] 



du 
u 2 - 1 



F(z,u), (95) 
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with the auxiliary function 
2 + 



F(z,u) 



sVT 

A 

H — r 
u 



-A 



~2uz 



2u 



3s 2 



(96) 



( Ae -2uz 



sVi + f 

l) coth [z — ln7 c (s)] . 

(97) 

We display in the main plot of Fig.[6]the potential profile 
Eq. (95 1 for various values of s. One sees in this plot that 



the behavior of the potential is mainly characterized by 
an interpolation between two regimes where the function 
tpu(z) changes its sign. In the first regime of large s (or 
weak surface charges and high salt concentrations) char- 
acterized by a positive energy barrier Sv (z) (see Fig. [5]), 
the resulting interfacial ion depletion responsible for a 
local ionic screening deficiency of the external potential 
increases the amplitude of the negative mean-field poten- 
tial <p(z). This regime is the reminiscent of the behavior 
observed in the previous part for the electrolyte system 
in contact with a dielectrically discontinuous wall (see 
Fig. [4]^a,) ) . Loosely speaking, in this regime, the MF the- 
ory overestimates the screening of the external potential. 
In the second regime of strong surface charges or dilute 
electrolytes corresponding to small values of s and an at- 
tractive self energy (see again Fig. [5]), the compact coun- 
terion layer formation at the interface is associated with 
a local ionic screening excess with respect to the bulk so- 
lution. This means that the MF theory underestimates 
the ionic screening in this parameter range, which is cor- 
rected with a positive tpu(z). The interpolation between 
both regimes will be quantitatively studied below. 

Expanding the function F(z,u) in Eq. (96 1 in inverse 



powers of s and carrying out the integral in Eq. ( 95 ) at 
z = 0, the surface potential follows in the form of a series 
as 

a 2 

tp u (0) = *L-r {-0.227s" 1 + 0.856s" 3 - 0.645s" 5 (98) 



-0.539s" ' -0.473s 



9 + O (s~ n 



)} 



From the relation (98), one finds that the one-loop cor- 
rection to the surface potential vanishes at the value 
s ~ 1.75 and remains negative for higher values of s. 
This is illustrated in the inset and the main plot of Figj6] 
The function ipu(0) also exhibits a minimum located at 
~ 3.19 before decaying to zero with increasing s. 

For large distances from the charged interface z> 1, 



the one-loop correction to the external potential Eq. ( 95 ) 
takes the simple asymptotic form 
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FIG. 7: (Color online) One-loop correction to counterion 
(solid lines) and coion densities (dashed lines) from Eq. ( 102 1 



for (a) s = 1000, (b) s = 1.5, and (c) s = 0.75. The inset 
displays the ionic self energy and the one-loop correction to 
the external potential for the same parameters. 



^i(«)~|-r 7c ( S )I( S )e-, 



where we introduced the auxiliary function 
I(s) = 



°° du 

2 - 1 



a" 



2 + s z 



- i 



-A 



2u 



3s 2 



(99) 



(100) 



The asymptotic form in Eq. (99 1 is displayed in Fig. [6] by 
black dashed lines. The function behind the exponential 
in this equation can be evaluated for large values of s as 



lc(s)l(s) 



-0.153s" 1 + 0.743s" 

„-7 i r\ 



0.628s" 5 



-0.542s"' +0(s" y ) . (101) 
One finds that this function vanishes at the particular 
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value s ~ 2, that is, the one-loop correction to the exter- 
nal potential changes its sign and becomes overall pos- 
itive at p, = 2/-t fc ~ 1 , before reaching a minimum located 
at s ~ 3.62. We note that interestingly, this sign re- 
versal takes place at the characteristic surface charge 
u s = \J pbj (27t£b), which is twice lower than the sur- 
face charge where the self energy Svq(z) becomes totally 
attractive. 

The deviations of the one-loop density from the MF 



density p± f (z) = p^e 



can be expressed as 



Ap±(z) 



pW(z) =-^ S M*)tMz). 



(102) 

We illustrate in the main plots of Fig. [7] the one-loop 
correction to ion densities renormalized by the coupling 
parameter T for s = 0.75, 1.5, and 1000. The one-loop 
corrections Svq(z) and if>u(z) are also shown in the inset 
of the same figures. In the case s — 1000 corresponding to 
weak surface charges where the amplitude of the poten- 
tial ipu(X) remains vanishingly small with respect to the 
repulsive self energy Svq(z) (see the inset of Fig.^a)), the 
repulsive solvation force induced by the interfacial salt 
screening loss is the only effect in play. Consequently, 
the one-loop correction lowers the MF density of both 
coions and counterions. 

In the case s = 1.5 of Fig. ^b) corresponding to a 
stronger surface charge or lower salt density where the 
interfacial counterion layer becomes dense enough to give 
rise to a strongly attractive branch of the self energy po- 
tential 5vq(z), correlation effects increase the MF density 
of both types of ions. In this range of the parameter s 
where the potential if>u{z) has a positive and a consid- 
erably large amplitude, the main plot of Fig. ^b) shows 
that interestingly, the additional attraction induced by 
the ionic self energy is amplified by the external poten- 
tial correction for coions, whereas the same attraction is 
partially cancelled for counterions. We show in Fig. 
that this effect is even stronger in the range s < 1 where 
the amplitude of the potential ipn(z) becomes compara- 
ble with the amplitude of the self energy 8vq(z). This 
indicates that correlations attenuate in this parameter 
regime charge separation in the interfacial region. 

We also see in Fig. [7|c) that for s < 1, the contribu- 
tion from the one-loop correction to the external poten- 
tial leads to a reduction of the MF counterion density at 
large distances from the interface. This peculiarity can be 
easily understood by comparing Eq. ( 99 ) with the large 
z asymptotic limit of the ionic self energy Eq. ( 94 1 . One 



sees that the one-loop correction to the external poten- 
tial ipu (z) is longer ranged than the self energy potential 
Svo(z). Thus, in the asymptotic limit i> 1, the former 
brings the main contribution to the one-loop corrections 
for the MF density in Fig.^c). 

Taking now into account the large distance asymptotic 
limit z > 1 of the MF potential Eq. ([38]) that can be 
written as ip(z) — —4j c (s)e~ z , the total one- loop exter- 
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FIG. 8: (Color online) Charge renormalization factor against 
s _1 for several values of the coupling parameter F. 



nal potential 4>u{ z ) = fi 2 ) + ^Pu{ z ) reads for z>l 



(103) 



where we introduced a charge renormalization factor 
dressed by electrostatic correlation in the form 



ri(s) = 2s 7c (s) 



1 -?K») 



(104) 



We illus trate in Fig. [8] the charge renormalization fac- 



tor Eq. (104) against the parameter s oc a s for differ- 



ent values of the coupling parameter T. The first point 
to be noted in this plot is the intersection between the 
curves for different values of T at s = 2 (see the inset). 
This point corresponds to the parameter range where the 
asymptotic large z limit of the one-loop correction to the 
external potential vanishes and all curves collapse onto 
the MF charge renormalization factor. One also sees that 
for s > 2 (the left portion of the intersection point), rj{s) 
increases with T. This behavior can be easily understood 
by noting that this parameter regime was shown above 
to correspond to a weak interfacial ionic screening defi- 
ciency resulting in a negative one-loop correction to the 
negative mean-field potential. Furthermore, by taking 



the vanishing surface charge limit of Eq. ( 104 1 , one finds 
linis^oo r)(s) — 1+3.8 x 10 -2 q 2 T. In other words, electro- 
static correlations interestingly yield a finite correction to 
the charge renormalization factor even in the limit of a 
vanishing fixed charge distribution. 

In the second regime s < 2 corresponding to stronger 
surface charges, it is seen that rj(s) changes its trend and 
starts to decrease with increasing coupling parameter F. 
Indeed, it was shown above that this regime is charac- 
terized by a compact counterion layer associated with 
an interfacial screening excess and a positive one-loop 
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correction to the negative MF potential. This aspect ex- 
plains the trend of the charge renormalization factor for 
s < 2. 

Furthermore, one notices that at a particular value of 
s, the factor ?y(s) changes its sign and becomes negative, 
resulting in a si gn reversal of the total one- loop poten- 
tial in Eq. ( |103[ ), which becomes positive. The reversal 
of the sign of the external potential above a character- 
istic surface charge is a signature of the charge inver- 
sion phenomenon. We note that this effect was also 
observed in Ref. [3D] for the symmetrically distributed 
electrolyte system. However, we also emphasize that the 
coupling parameter range where this effect takes place in 
Fig. [8] is beyond the validity of the one-loop approxima- 
tion. Indeed, we verified that in this electrostatic cou- 
pling regime, even the perturbative solution scheme of 
the SC equations ([6|-([7]) introduced in section III A does 
not converge. Furthermore, the MF counterion densi- 
ties reach in this charge density regime unrealistic values 
beyond the close packing, an artefact known to origi- 
nate from the absence of ionic excluded volume effects in 
the model Hamiltonian of Eq. (JlJ . The inclusion of 
hard-core effects known to reduce the interfacial counte- 
rion densities is also expected to shift the charge reversal 
point to larger surface charges. Consequently, it is clear 
to us that in Fig. [HJ the curves with T > 1 have no quan- 
titative reliability for large values of s _1 . That being 
said, we believe that this result still presents some quali- 
tative interest since the effect observed in Fig.[8]might be 
the precursor of the actual charge inversion phenomenon 
in asymmetrically distributed salt systems. This point 
needs however to be verified by comparisons with MC 
simulations. 



V. SUMMARY AND CONCLUSION 

In this article, we investigated electrostatic correlation 
effects in symmetric electrolytes in contact with charged 
planar interfaces separating the solvent region from a 
membrane area free of ions. To this aim, we introduced in 
the first part of the work two computational approaches 
to solve the electrostatic SC equations derived in Ref. [TS] 
in the presence of dielectric discontinuities where the one- 
loop theory fails. Then, we compared in the second part 
the theoretical ion density profiles obtained from these 
computation schemes with the results of our MC simu- 
lation data in order to determine the validity domain of 
the SC equations at neutral dielectric interfaces. It was 
shown that the DH theory that neglects the interfacial 
variations of the ionic screening exhibits a quantitative 
accuracy up to the characteristic bulk density pt ~ 0.01 
M, while the SC theory remains accurate up to pi, ~ 0.2 
M, thus improving the quantitative accuracy of the DH 
theory by one order of magnitude in ionic strength. The 
deviations of the SC results from MC simulation data at 
Pb 0.2 M may be either due to electrostatic correlation 
effects, or hard-core effects that become relevant in this 



concentration regime. This point can be enlightened in 
a future work by solving the extended SC equations of 
Ref. 23 that account for the excluded volume of ions. 

Within the validity regime of SC equations, we also 
validated the accuracy of a restricted variational scheme 
introduced in [5T] for slit nanopore systems. This obser- 
vation may explain the success of similar methods fre- 
quently used in nanofiltration studies in order to predict 
experimental salt rejection rates [IOlU2j . Furthermore, 
we computed within the WKB formalism the interaction 
energy between a charged rigid polymer and a neutral di- 
electric interface. We showed that due to the interfacial 
ionic screening deficiency neglected in the DH approxi- 
mation, the energetic cost to bring a polymer from the 
bulk to the proximity of the interface is characterized 
by a different scaling law with the distance from the in- 
terface as well as a significantly higher amplitude than 
the prediction of the DH theory j^Sj. This observation 
relevant for protein-surface interactions could be verified 
with MC simulations. 

In the third part of the article, we considered the case 
of ions at charged interfaces. For ions in the proximity of 
a weakly charged dielectric wall, we showed that the main 
correlation effect is a strong interfacial salt exclusion re- 
sulting in an ionic screening deficiency that increases the 
amplitude of the negative MF external potential. The 
latter effect strengthens the counterion attraction and 
coion repulsion far away from the interface. Thus, corre- 
lation effects in weakly charged membrane nanopores are 
expected to amplify the charge separation phenomenon 
induced by Donnan exclusion. Comparing the ion den- 
sity and external potential profiles obtained from the SC 
scheme with MC simulation data of Ref. [51], we also 
showed that the SC theory is able to accurately handle 
these correlation effects in a parameter regime where the 
MF theory exhibits a significant deviation from the MC 
results. 

The forth part of the work dealt with electrostatic cor- 
relation effects in dielectrically homogeneous systems at 
the one-loop level. We first expanded the SC theory at 
one-loop order and found that one-loop corrections in this 
system are driven by the competition between the inter- 
facial salt screening loss driving the ions towards the bulk 
area, and the counterion screening excess induced by the 
interfacial counterion layer attracting them to the inter- 
face. This competition can be quantified in terms of the 
balance between the characteristic thickness of the inter- 
facial counterion layer p and the radius of the ionic cloud 
k^ 1 around a central charge in the bulk region. Namely, 
in the presence of a weak surface charge corresponding to 
a diffuse counterion layer p ^> k^ 1 , the surface-repulsive 
salt screening effect is the dominant mechanism. In this 
case, correlation effects result in a net interfacial exclu- 
sion of both coions and counterions, and a weak negative 
correction to the MF potential. 

Decreasing the bulk salt concentration or equiva- 
lently increasing the surface charge to the range a s ~ 
yj2pb/ (ttIb), the thickness of the counterion layer p be- 
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comes comparable with the ionic cloud radius n^ 1 . As a 
result, the surface-attractive counterion screening effect 
starts to set on and correlation effects increase in this 
regime the MF density of both types of ions. Moreover, 
due to the strong counterion excess close to the surface, 
the one-loop external potential correction also acquires 
in this surface charge regime a positive and large am- 
plitude. Consequently, while a further increase of the 
surface charge will result in a amplification of the MF 
coion density in the whole half space z > and the MF 
counterion density close to the surface, the latter will be 
attenuated by correlation effects outside the interfacial 
region. 

It was also shown that if one reaches a high enough 
surface charge value, the one-loop theory predicts charge 
inversion for T > 1. However, as we stressed in the main 
text, this result should be considered with caution. In- 
deed, this parameter regime stays well beyond the valid- 
ity range of the one-loop theory. Hence, comparison with 
MC simulations is needed in order to verify whether this 
observation might be the precursor of the actual charge 
inversion effect in asymmetrically distributed salt sys- 
tems. 

The one-loop theory of asymmetrically partitioned 
electrolytes presented in the final part of this work 
bridges a gap between the DH theory of symmetric 
salts at neutral interfaces [55] and the one-loop theory 
of counterions in contact with a charged surface |15j . 
Furthermore, we note that the theoretical concepts 
introduced in the present work can be applied to 
more complicated geometries such as cylindrical ion 
channels [52], or the charged Yukawa model studied 
in Ref. [23] in order to evaluate the importance of the 
ionic excluded volume on the correlation effects that we 
have investigated. We would also like to establish in a 
future work the validity regime of the SC equations Q 
and ([7]) with respect to the surface charge by running 
extensive MC simulations of ions in contact with charged 
interfaces. Finally, we note that the perturbative SC 
scheme can also account for electrostatic correlation 
effects in more complicated electrostatic systems where 
the one-loop theory fails, such as polymer brushes at 
charged dielectric interfaces or dielectric polyelectrolytes 
in ionic solutions. We would like to treat these cases in 
future works. 



Appendix A: Perturbative solution of the closure 
equations within WKB approach 



A perturbative solution to the closure equations (54) 



and ( 66 1 was derived in Ref. [5J for simple interfaces in 
the case e m = 0. We will rederive in this Appendix this 
result obtained for finite e rn . The first step consists in 
writing the screening function of Eq. ( 54 ) in the form 



[1 + Sn (z)} 



(Al) 



where Sno(z) is a "small" perturbation of the screening 
length induced by the image potential that we choose as 
the screened image potential within the undistorted ionic 
atmosphere approximation |10j . 



dno(z) = exp 



q 2 e B A 2Kb 

Az 



1. 



(A2) 



We note that this choice is motivated by the very weak 
dependence of the image potential profile on the input 
image potential, a good convergence property of the it- 
erative computation scheme observed in the numerical 



solutions of the closure equations ( 54 ) and ( 66 ) . Further 



more, the naught in Eq. (|AT]) denotes the input function 
at the zeroth order iterative level. By injecting this func- 
tion into the electrostatic potential (67) and expanding 
in Sn Q (z), one gets after some algebra 



6v(z) = £j ^e~ 2KbZ +Svi(z), 
with the correction to the WC image potential 



5vi{z) 



e B Kb 

Ib 



1 



exp 



-2Kb -S 



(A3) 



(A4) 



K b A 5n (z)e 



+£ B A K 2 b T(Q,2 Kb z) 



a j zl d5n (z') 
dz' 



where Y(a,x) is the incomplete Gamma function |35) . 
By injecting now Eq. (A3) into Eq. the ion density 
profile at the first iterative level finally takes the form 



Pb 



exp 



-2KbZ 



4z 



Svi_(z) 



(A5) 
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Appendix B: Third order perturbative solution of 
SC equations for neutral interfaces 

We will give in this Appendix the third order pertur- 
bative correction for the solution of the SC equation ^ 
for a vanishing surface charge a(z) = 0. Because this 
case corresponds to a zero external potential, the com- 
ponents of the correction to the Green's function Sv nm 
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with m > vanish. As in the computation of the second 



Appendix C: Monte Carlo simulations 



order calculation presented in section III A the analytical 
task consists in injecting into Eq. (17) the expansion of 



the Green's function ( 19 ) and keeping only the terms up 
to X 3 ,. One obtains, in addition to the first and second 



order corrective terms in Eqs. (21) and (23), the third 
order correction in the form 



Sv 30 {z) = -~^Pb I dz 1 n (zi) ^Svl^zx) - 45u 20 (zi)] 



-V w (z 1 ) 



x/ 2 (r,r', zi)e 



-2q p b J dzidz2n (zi)Svio(zi)Sn a (z 2 ) 

x [h(r,r',z 1 ,z 2 ) + h(r' ,t,z u z 2 )\ 

Xe -V w {z 2)] 

-8q 6 pl J dzidz 2 dz 3 8no(zi)Sno(z 2 )Sno(z 3 ) 
x/ 4 (r,r',z 1 ,Z2,z 3 ), (Bl) 
where we defined the new function 

h(v,r',z 1 ,Z2,z 3 ) = J ^e' k -( r ir r i)i„( ZlZl ^) 

xv (zi,Z2,k)v (z 2 ,Z3,k) 
xv (z 3 ,z',k). (B2) 

Finally, at the third order perturbative level, the ion den- 
sity Eq. ^ reads 



p(z)/p(°Hz) = l-\ v ^6v 10 (z) 

a 2 
,2<7 



(B3) 



[Svf (z)~4Sv 20 (z)] 

"A^ [q 4 Sv 3 w (z) - l2q 2 Sv w (z)Sv 20 (z) 
+246v 30 (z)} . 



We present in this Appendix the details of the canon- 
ical MC simulations. During the simulations, the total 
particle number was fixed at N p ~ 4096, with the cor- 
responding particle density pmg = N p /(l x l y l z ) = p b , the 
size of the simulation box l x — l y — (Np/l.bpb) 1 ^ 3 and 
l z = 1.5Z X , and impenetrable walls located at z = and 
z = l z . Periodic boundary conditions in the x and y di- 
rections were used. The value of l z that we considered 
was large enough so that the box geometry was equivalent 
to the actual single interface geometry in consideration. 



The MC simulations were run on nvidia graphical pro- 
cessing units (GPU) using CUDA algorithm [39]. In the 
simulations, all the relevant variables were set to reside 
in GPU RAM memory. By dividing the work in small 
independent parts that are all executed in parallel, the 
GPU can provide up to two orders of magnitude improve- 
ment in performance. In our system, we have divided the 
calculation of the energy difference between the present 
and the trial configuration in small tasks. For each trial 
move, 27Vp threads were used for calculating the indi- 
vidual particle-particle interactions, followed by a sum- 
reduction to obtain the total energy difference between 
the present and the new configuration. Finally, only 1 
thread was checking the acceptance conditions and mak- 
ing the necessary updates. Even if some parts of the MC 
algorithm were not able to fully utilize the GPU, we were 
able to obtain a speed-up larger than 100 times. 
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